set tcl_precision 8

puts " "
puts "======================================================="
puts "=          Interaccion magnetica Confinado            ="
puts "======================================================="
puts " "
puts "[code_info]"
set pi 3.14159265359

#############################################################
# ENTRADAS DE LA SIMULACION				    #
#############################################################
#Reiniciar simulacion    				    
set restart "yes"
#Flujo activado?
set fluxon "no"
#Analisis						    
set setAnalisis "no"
#Sistema confinado					    
set setConfined "yes"
#Imprimimos las configuraciones				    
set printConfig "yes"
# Interaccion magnetica activada?
set magon "yes"

#Paso de tiempo 					    -
set dt 0.0005
set dth [expr 0.5*$dt]
setmd time_step $dt
#Coeficiente de friccion de ec. Langevin
set gamma 10.0
#Temperatura
set temp 1.0
#Ratio lbox_z/lbox_x (Se considera lbox_x=lbox_y)	    
set ratio 0.5
#Fraccion en volumen					    
set phi 0.05
#Intensidad del campo magnetico				    
set u0 [expr 100.0/8.0]
set m [expr sqrt($u0)]
puts "U0: $u0 (m: $m)"

#Intervalo de impresion de magnitudes energeticas
set stepE 1000

if {$fluxon == "yes"} {
    #solamente necesario si se activa el flujo
    #Amplitud de la deformacion senoidal en tanto por uno
    set gamma_cero 0.002
    #Omega
    set omega 10.0
    #Duraccion de la simulacion y frecuencia de almacenado
    #Numero de configuraciones por ciclo almacenadas
    set nconfciclo 100
    #Numero total de ciclos almacenados
    set nciclototal 100
    #Paso final de la simulacion
    set naccfinal [expr int(100*(2.0*$pi)/($omega*$dt))]
    set ialmacen [expr int(2.0*$pi/($omega*$dt*$nconfciclo))]
    #Frecuencia de almacenamiento de las configuraciones
    set nsave $ialmacen
}

set nFileRestart "system_ini.dat"
#Numero de particulas					    
set npart  0

if { $restart == "no" } {
    set volume 0.0
    set objFile [open "posicion.dat" "r"]
    set datos [read $objFile]
    close $objFile
    set records [split $datos "\n"]
    set linea(0) 0
    set maxdiam 0.0
    
    ## Iterate over the records
    foreach rec $records {
	    ## Split into fields on colons
	    set fields [split $rec " "]
	    set tonto ""
	    set n 0
    
	    foreach pp $fields {
		    if { $pp != ""} {
			    set linea($n) $pp
			    #puts "$n $pp"
			    incr n
		    }
	    }
	    set x $linea(0)
	    set y $linea(1)
	    set z $linea(2)
	    
	    set diam $linea(3)
	    if {$diam > $maxdiam } {
		    set maxdiam $diam
	    }
	    set volume [expr $volume + $diam*$diam*$diam*$pi/6.0]
	    incr npart
    }
	puts $volume
    set npart [expr $npart -1]
	    
}

if { $restart == "yes"} {
	set volume 0.0
	set objFile [open "system_ini.dat" "r"]
	set datos [read $objFile]
	close $objFile
	set records [split $datos "\n"]
	set linea(0) 0
	set maxdiam 0.0
	
	## Iterate over the records
	foreach rec $records {
		## Split into fields on colons
		set fields [split $rec " "]
		set tonto ""
		set n 0
	
		foreach pp $fields {
			if { $pp != ""} {
				set linea($n) $pp
				#puts "$n $pp"
				incr n
			}
		}
		set x $linea(0)
		set y $linea(1)
		set z $linea(2)
		
		set vx $linea(3)
		set vy $linea(4)
		set vz $linea(5)
		
		set fx $linea(6)
		set fy $linea(7)
		set fz $linea(8)
		
		set diam $linea(9)
		if {$diam > $maxdiam } {
			set maxdiam $diam
		}
		set volume [expr $volume + $diam*$diam*$diam*$pi/6.0]
		incr npart
	}
	set npart [expr $npart -1]
}

#Volumen de la caja					  
set volume     [expr $volume/$phi]
#Lado de la caja					    -
set lx  [expr pow($volume/$ratio,1.0/3.0)]
set ly $lx
set lz [expr $ratio*$lx]
puts "Lado de caja, lx: $lx $ly $lz"

#############################################################
# INICIALIZACION                                            #
#############################################################

if {$setConfined == "yes"} {
	#Confinamiento en eje z					    -
	puts "Hay confinamiento (barreras en eje z)"
	setmd periodic 1 1 0
} else {
	#No Confinamiento 					    -
	puts "No hay confinamiento"
	setmd periodic 1 1 1
}

#Creamos la caja 					    
setmd box_l $lx $ly $lz
#Numero de ciclos 	
if {$fluxon == "yes"} {
	set n_cycle $naccfinal
} else {
	set n_cycle 100000000
}

#Numero de pasos por ciclo				    
set n_steps 1
#Para la lista de vecinos 				    
setmd skin 0.5
#Simulacion NVT			                            
integrate set nvt		                            
thermostat langevin $temp $gamma
#Archivo de estado de la simulacion                         
set nFileSystem "system_bk.dat"
#Impresion de resumen de la simulacion                      
puts [setmd box_l]
puts [setmd time_step]
puts [integrate]
puts [thermostat]

#############################################################
# PREPARAMOS LOS CORRELADORES                               #
#############################################################

if {$setAnalisis =="yes"} {
    set base 1.25
    set npc 32
    set nc [expr 1000000/$npc]
    set k 0
    set kp 0
    set tciclo 0
    
    for {set i 0} {$i<500000} {incr i} {
    	set tempi($i) 0
    	set nss($i) 0
    }
    
    for {set i 0} {$i < $nc} {incr i} {
    	for {set j 0} {$j <= $npc+1} {incr j} {
    		set val [expr $tciclo + pow($base, $j)]
    		set tempi($k) val
    		set nss($kp) [expr int($val)]
    		if {$kp >=1 && $nss($kp) == $nss([expr $kp -1])} {
    			set kp [expr $kp -1]
    		}
    		incr k
    		incr kp
    	}
    	set tciclo $nss([expr $kp -1])
    }
    
    set inext $nss(0)
    set icorr 0
	
}

#############################################################
# INTERACCIONES                                             #
#############################################################
#set cut 2.7
#set sig 0.1
#set r0 0.9
#set cut [expr $r0+2.5*$sig] 
#inter 0 0 lennard-jones 1.0 $sig $cut auto $r0
#inter 0 1 lennard-jones 1.0 0.5 [expr pow(2.0,1./6.)] auto 0.0

# Para interaccionees sin atracciones con potencial klingemberg
set cut 4.0
inter 0 0 klingemberg $u0 100.0 $cut
inter 0 1 klingemberg $u0 100.0 [expr $cut]


#############################################################
# CREAMOS LAS PARTICULAS                                    #
#############################################################

if {$restart == "no"} {
    set objFile [open "posicion.dat" "r"]
    set datos [read $objFile]
    close $objFile
    set records [split $datos "\n"]
    set linea(0) 0
    set nnn 0
    ## Iterate over the records
    foreach rec $records {
	## Split into fields on colons    
	set fields [split $rec " "]
	set tonto ""
	set n 0
	foreach pp $fields {
	    if { $pp != ""} {
		    set linea($n) $pp
		    incr n
	    }
	
	    set x $linea(0)
	    set y $linea(1)
	    set z $linea(2)
	    set diam $linea(3)
	    set masa [expr $diam*$diam*$diam]
	    if {$nnn >= $npart} {
		    set n $npart
	    }
	    part $nnn pos $x $y $z dipm 1.0 dip 0 0 [expr $m*$masa] diameter $diam mass $masa type 0
	}   
	incr nnn
    }
    set min [analyze mindist]
    puts "Distancia minima: $min"
}
	
if { $restart == "yes"} {
    set volume 0.0
    set objFile [open "system_ini.dat" "r"]
    set datos [read $objFile]
    close $objFile
    set records [split $datos "\n"]
    set linea(0) 0
    set nnn 0
    
    ## Iterate over the records
    foreach rec $records {
	## Split into fields on colons
	set fields [split $rec " "]
	set tonto ""
	set n 0
	
	foreach pp $fields {
	    if { $pp != ""} {
		set linea($n) $pp
		#puts "$n $pp"
		    incr n
	    }
	}
	set x $linea(0)
	set y $linea(1)
	set z $linea(2)
	
	set vx $linea(3)
	set vy $linea(4)
	set vz $linea(5)
	
	set fx $linea(6)
	set fy $linea(7)
	set fz $linea(8)
	
	set diam $linea(9)
	set masa [expr $diam*$diam*$diam]
	if {$nnn < $npart} { 
	    part $nnn pos $x $y $z v $vx $vy $vz f $fx $fy $fz dipm 1.0 dip 0 0 [expr $m*$masa] diameter $diam mass $masa type 0
	}
	incr nnn
    }
    set min [analyze mindist]
    puts "Distancia minima: $min $npart"
}

#############################################################
# CONFINAMIENTO                                             #
#############################################################
if {$setConfined == "yes"} {
    #Creamos dos paredes				    
    constraint wall normal 0 0 -1 dist   [expr -1.0*$lz] type 1
    constraint wall normal 0 0 1 dist  [expr 0.0] type 1
}

#############################################################
# INTERACCION MAGNETICA                                     #
#############################################################
#Interaccion magnetica vecinos cercanos  
if {$magon == "yes"} {
	inter magnetic 1.0 dawaanr
}

# CREAMOS ARCHIVOS DE OBSERVABLES	
set obsEner [open "energia.dat" "w"]

#############################################################
# BUCLE PRINCIPAL			                    #
#############################################################

set i 0
set ncount 0
set idconfig 0
set k 0
while { $k<$n_cycle } {
    set ncount 0
    set t [expr $k*$n_steps*$dt]
    if {$fluxon == "yes"} {
        couette $gamma_cero $omega $t 1.0
        
        #Tiempo al que corresponde la configuracion
        if {$k%$nsave == 0} {
        	printConfig
        }
    }
    integrate $n_steps
    # ANALISIS 			                    
    if {$setAnalisis =="yes"} {
        if { $k >= $inext } {
		puts "TIEMPO: $t"
		printConfig
        	incr icorr
        	set inext $nss($icorr)
        }
    }
	
    # IMPRIMIMOS COSAS		   
    if {$k%$stepE == 0} {
        set emag [analyze energy magnetic]
        set ekin [analyze energy kinetic]
        set etot [analyze energy total]
        set erep [analyze energy nonbonded 0 0]
        #set erep [expr $etot-$emag-$ekin]
        set temp [expr $ekin/((3.0/2.0)*$npart)]
        set hora [clock format [clock seconds] -format "%H:%M:%S"]
        puts "($hora) Time: $t, Etotal: [expr $etot/$npart], Temp: $temp"
	if {$magon == "yes"} {
		puts $obsEner "$t [expr $emag/$npart] [expr $erep/$npart] [expr $ekin/$npart] [expr $etot/$npart]"
	} else {
		puts $obsEner "$t [expr $erep/$npart]  [expr $ekin/$npart] [expr $etot/$npart]  "
	}
        flush $obsEner
        #puts "Emagn: $emag"
        set systembk [open $nFileSystem "w"]
        puts $systembk "treal $t u0 $u0  npart $npart lbox $lx $ly $lz"
        if {$magon == "yes"} {
	    puts $systembk "[expr $emag/$npart] [expr $erep/$npart] [expr $ekin/$npart] [expr $etot/$npart] "
        } else {
	    puts $systembk "[expr $erep/$npart] [expr $ekin/$npart] [expr $etot/$npart] "
        }
        for {set i 0} { $i < $npart } {incr i} {
        	      puts $systembk [part $i print pos v f diameter]
        }
        close $systembk
    }
    
	incr k
}

close $obsEner
